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ABSTRACT 

We address the question of whether the formation of high-mass stars is similar to 
or differs from that of solar-mass stars through new molecular line observations and 
modeling of the accretion flow around the massive protostar IRAS20126+4104. We 
combine new observations of NH3 (1,1) and (2,2) made at the Very Large Array, new 
observations of CH3CN(13-12) made at the Submillimcter Array, previous VLA obser- 
vations of NH3(3,3), NR"3(4,4), and previous Plateau de Bure observations of C 34 S(2- 
1), C 34 S(5-4), and CH 3 CN(12-11) to obtain a data set of molecular lines covering 
15 to 419 K in excitation energy. We compare these observations against simulated 
molecular line spectra predicted from a model for high-mass star formation based on 
a scaled-up version of the standard disk-envelope paradigm developed for accretion 
flows around low-mass stars. We find that in accord with the standard paradigm, the 
observations require both a warm, dense, rapidly-rotating disk and a cold, diffuse in- 
falling envelope. This study suggests that accretion processes around 10 M stars are 
similar to those of solar mass stars. 

Key words: Keywords. 



1 INTRODUCTION 

Does the formation of massive stars differ significantly from 
that of solar mass stars? As far as we know, stars of all 
masses form in gravitationally unstable regions of molec- 
ular clouds and gain their mass by accretion. A standard 
model developed for accretion flows around low-mass stars 
consists of two-components, a rotationally-supported disk 
inside a freely-fa lling envelope (|Shu. Adams. Lizand 1 19871 ; 
lHartmanrJ 2001). This model has been particularly suc- 



cessful in explaining infrared observations of low-mass star 
formation. The disk and envelope produce an excess of 
long-wavelength infrared emission that has been adopted as 
the identifying signature of accreting proto stars in Galac- 
tic (| Whitney et al.ll2003l ) and extragalactic l|Whitnev et al.l 
120081 ) star-forming regions. Furthermore, because the disk 
and envelope have different densities and temperatures, the 
evolutionary state of the protostars can be identified by the 
shape of the infrared spectral energy distribution: c lass 
(envelope dominated) a nd class I (disk dominated) (|Ladal 
ll987l ; lAndre et al.lll993T ). 

Does this standard two component accretion model de- 
veloped for low-mass stars also describe the accretion flows 
around massive stars? There are some doubts. The more 



massive stars are luminous enough to generate radiation 
pressure and hot enough to ionize their own accretion flows 
such that the outward radiative and th ermal pressures rival 
the inward pull of the stellar gravity llLarson fc Starrfleldl 
Il97ll ; lKahnlll974l ; lKetdl2002l ; iKeto fc Woodll200d ). Do these 
outward pressures result in accretion flows that are different 
around more massive stars? In this paper we compare new 
and previous molecular line observations of an accretion flow 
around one massive star against the standard disk-envelope 
paradigm for accretion flows around low-mass stars. 

The previous molecular line observations of the mas- 
sive protostar IRAS20126+4104EI, suggest an accretion disk 
and bipolar outflow around a 7 to 1 5 protostar em- 



bedded in a dense molecular e nvelope (ICesaroni et al.l 
Zhang. Hunter, fc SridhararJ 19981: Cesaroni et al.l 



1997 



1999 



iKawamura et al. 1999: Zhang ct al. 1999; Shepherd ct al 
hood: ICesaroni et al.l ]2005l; (Lebron et alj hood; ISu e~ 



EMM 

iu et al.ll2008l ). Previous infrared observations of ab- 
sorption and scattering also r eveal a disk and outflow cavity 
immediately around t he star l|Sridharan. Williams, fc Fuller! 
120051 : IdeBuizer |[2007l ) . 

In this study, we assemble a suite of observations of 
molecular lines of different excitation temperature in or- 
der to compare with molecular spectra predicted from the 
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1 IRAS20126 is located in the Cygnus-X region at a distance of 
1.7 kpc (Wilking et al. 1989). 
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disk-envelope model. Lines of different excitation temper- 
ature are useful because massive stars heat the surround- 
ing molecular gas to observationally significant temperatures 
(> WOK) at observationally significant distances from the 
star, (Tgas ~ T*(-R/i?*) a ), and we can exploit the relation- 
ship between temperature and radius to distinguish emis- 
sion from gas at different radii around the star. We expect to 
identify the emission from the higher excitation temperature 
lines with gas in the flow closer to the star and also to sep- 
arate the emission from t he disk and envelope compon ents. 
In previous ob s ervat ions, iKeto. Ho fc Haschickl (|l987f ) and 
ICesaroni et al.l (1 19941 ) used this technique, observing several 
lines of NH3 to study the accretion flows around very high 
mass stars associated with HII regions. 

We present new observations of the NH3 (1,1) and (2,2), 
inversion transitions made with the National Radio As- 
tronomy Observatory's Very Large Array (VLA)0 and new 
observations of CH3CN(13-12) made with the Submillime- 
ter Array (SMA). The new observations of NH3(1,1) and 
Nti3(2,2) have a factor of 3 better sensitivity than th e earlier 
observations of IZhang. Hunter, fc Sridharanl (1 19981 ) . Com- 
bined with previou s observations of NH3(3,3) and NH3(4,4) 
<|Zhang et al.lll999h the 4 NH3 lines span a range of excita- 
tion temperatures from 23 K to 200 K. We also have addi- 
tional lines of lower and higher excitation temperature from 
previous observations of the C 34 S(2-1) and C 34 S(5-4) lines 
with energies of 7 K and 35 K respectively, and the ladder 
of CH3CN (J=12-ll) lines with energies ranging from 69 K 
for (K=0) to 419 K for (K=7). The C 34 S and CH 3 CN obser - 
vati ons were previously p resented in lCesaroni et alj (|l999l ') 
and lCesaroni et all (|2005h . 

We specify the two-component accretion model in terms 
of 6 parameters: the scale factors for the density and tem- 
perature of the envelope and of the disk, the angular mo- 
mentum of the envelope, and the stellar mass. We use our 
molecular line emission code MOLLIE to predict molecular 
line spectra from the parameterized model, and we use a 
least-squares fitting procedure to adjust the parameters to 
fit the observations. 

We find that the two-component, disk-envelope model 
can successfully describe the observations, but a single com- 
ponent model cannot. A warm, dense rotationally-supported 
disk is required to obtain sufficient brightness and width in 
the high-excitation lines, and a cold, large-scale envelope is 
required to match the emission from the lower excitation 
lines. We find no evidence that the accretion flow around 
IRAS20126 is profoundly altered by the outward force of ra- 
diation pressure or by ionization. At 10 Mq, the star might 
simply not be luminous enough or hot enough for its ra- 
diation pressure or ionization to significantly affect the ac- 
cretion process. Based on this example, the accretion flows 
around 10 Mq stars are quite similar to flows around lower- 
mass stars. 



2 The National Radio Astronomy Observatory is a facility of the 
National Science Foundation operated under cooperative agree- 
ment by Associated Universities, Inc. 



2 THE DATA 

2.1 Ammonia inversion lines 

We observed the IRAS 20126+4104 region, R.A. (J2000) 
= 20:14:26.06, decl. (J2000) = 41:13:31.50, in the inversion 
transitions of NH 3 (J,K) = (1,1), (2,2), (3,3) and (4,4) on 
1999 March 27 and 1999 May 29 with the VLA in its D 
configuration. We used the 2IF correlator mode to sam- 
ple both the right and left polarizations, a spectral band- 
width of 3.13 MHz, and a channel width of 24 kHz or 0.3 
kms -1 . The primary beam of the VLA at the NH3 line fre- 
quencies is about 2'. Quasars 3C48 and 3C286 were used 
for flux calibration, 3C273 and 3C84 for bandpass cali- 
bration, and quasar 2013+370 for phase calibrations. The 
NH3(3,3) and NH3(4 ; 4) obs ervations were previously dis- 
cussed in IZhang et all (l999l ). 

The VLA data were processed using the NRAO As- 
tronomical Image Processing System (AIPS) package. Data 
from the two days were combined to achieve an rms of 3 
mjy per 0.3 km s -1 channel for the NH3 (1,1), (2,2) lines 
and 6 mjy and 4 mjy per 0.6 km s -1 channel for the NH3 
(3,3) and (4,4) lines respectively. Figure [T] shows two images 
of IRAS20126 in the integrated intensity of NH 3 (1,1) and 
(2,2). The protostar is located within the bright core of the 
larger scale cloud. Also shown in figure [T] is the image o f 
the bipolar outflow in SiO(2-l) from ICesaroni et al.l l|l999f ). 
The NH3 spectra used in the analysis are taken at positions 
offset from the phase center by -0.23, 1.00 (center), -3.10, 
3.25 (left), and 2.50,-0.94 (right) in arc seconds of RA and 
dec. The corresponding absolute positions are R.A. (J2000) 
= 20:14:26.02, decl. (J2000) = 41:13:32.8 (center) , R.A. 
(J2000) = 20:14:25.79, decl. (J2000) = 41:13:30.8 (right) , 
and R.A. (J2000) = 20:14:26.28, decl. (J2000) = 41:13:35.0 
(left) , The locations of the spectra are shown as red stars 
on figure [T] 



2.2 Methyl Cyanide and Carbon Sulfide lines 

The CH 3 CN(13-12) observations were made with the SMA 
on July 10, 2006 with a maximum baseline of ~ 0.5 km, a 
beam size of 0.40 x 0.36" , and a phase center of R.A. ( J2000) 
= 20:14:26.02, decl. (J2000) = 41:13:32.7. The channel width 
is 0.8 kms -1 at the observing frequency of 239 GHz. The 
rms noise of the data is 0.05 Jy beam -1 channel -1 or 8 K 
channel -1 in brightness temperature. 

Data in the CH3CN J = 12-11, and C 34 S J = 5-4 and 
J = 2-1 transitions were obtained from the IRAM interfer- 
ometer between 1996 December through 1997 March, and 
2002 January through March. The phase center of the ob- 
servations is R.A. (J2000) = 20:14:26.03, decl. (J2000) = 
41:13:32.7. The primary beam of the array at the frequency 
of the observed CH3CN transitions (230 GHz) and C 34 S(5- 
4) (241 GHz) is approximately 20" with a spatial resolu- 
tion of about 0.8". At the lower frequency of the C 34 S(2-1) 
transtion (96 GHz) the primary beam is about 50" and the 
spatial resolution is and 2.4". The full details of these ob- 
servations can be found in Cesaroni et al. (1999, 2005). 
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3 THE STANDARD MODEL FOR STAR 
FORMING ACCRETION FLOWS 

3.1 The density of the accretion flow 

The standard model for a star-forming accretion flow con- 
sists of an inflowing envelope around a ro tationally sup- 
orted disk. Similar to the approach in IWhitnev et al.l 
20031 ). we use the model of lUlrichl l| 19761 ) to describe the 
gas den sity in the inf alling envelope and standard thin-disk 
theory (|Pringldll98ll ) to describe the disk, 

The model for the envelope assumes that the gas flows 
toward a gravitational point source along ballistic trajecto- 
ries. The flow conserves angular momentum and ignores the 
pressure and se lf-gravity of the gas. T his is the same enve- 
lope model that lWhitnev et al.l (|2003l ) adopted for low-mass 
stars and similar to the model used to describe the accre- 
tion flow onto a cluster of high-mass stars, G10.6-0.07 jKetol 
Il99d : iKeto fc Woodll2006l ) . Thus despite, or perhaps because 
of, its simplicity, the model has found application to accre- 
tion flows from low-mass stars to star clusters. The model is 
particularly useful in analyzing observational data because 
only 3 parameters are required to describe individual cases: 
1) the density of the gas at a single radius (the gas density 
elsewhere follows from mass conservation), 2) the mass of 
the point source, and 3) the specific angular moment um. 

The model for the envelope is fully described in 
(|l976h ^ We use the equatio n s as p resented in iKetol 
and Mendoza. Canto. Ragal (|2004l ). The envelope density is 
(equation 10 of lKetol (|2007l ) ) , 



Pcnv = peo(r/Rd) 



-3/2 



cos 6 



-1/2 



x {I + (r / R d y 1 (3 cos 2 Oo-l))- 1 



(1) 



where p e o = Penv(Rd) is the density in the mid-plane, 
6 — 7r/2, at radius, Rd, #o is the initial polar angle of the 
streamline, and 8 and r are the polar angle and radius at 
each point along a streamline. The angle, 9 is related to the 
polar radius, r, and the initial angle do (equation 7 of lKetol 

1153)), 



Rd cos #o sin 0q 
cos 8o — cos 6 



(2) 



The density is re lated to the mass accretion rate (equation 
11 of lKetol l|2007l )) 



M = p e0 {Rd)^R d vk 



(3) 



where Vk is the Keplerian velocity at Rd. 

The model for the disk assumes that a rotationally- 
supported disk forms at the radius where the centrifugal 
force in the rotating envelope equals the gravitational force 
of the point mass, 

El = 9M* (4) 

R-l Rl () 



and F = v^R. The disk is truncated at Rd- 

The gas density in the disk is (equation 3.14 of IPringlc 



p disk (z,R) = p a (R)e X p(-z 2 /2H 2 



(5) 



where po(R) is the density in the mid-plane, 6 = 7r/2, at any 
radius R. We use R — \/ x 2 + y 2 to denote the cylindrical 



radius, and r = ^/ x 2 + y' 2 + z 2 for the polar radius. In the 
thin disk theory, the scale height, H, is 

H 2 = c 2 s R 3 /GM, (6) 

where c s is the sound speed. We follow lWhitnev et all ([20031 ) 
and use a modified version of this equation, 



H = H {R/R* 



(7) 



with Ho = 0.01i?». The density in the mid-plane of the 
accretion disk, po(R), is 



Po = Pdo(Rd/RY 



(8) 



where pdo is the densi t y in t he mid-plane at Rd. Again fol- 
lowing |whitnex!elaD (|2003l ), we use an exponent of 2.25 in 
equation [8] rather than 2 which is der i ved in the thin-disk 
theory. As explained in lWhitnev et afl (120031 ) the modifica- 
tions to the equations for the scale height and mid-plane 
densities are based on fits to numerical models of disk struc- 
ture. 

The densities in the disk and envelope at radius Rd are 
related by a factor, A p , 



PdO = A p p e o 



(9) 



For example, in the steady-state flow, the disk would have a 
higher density if, as should be the case, the inward velocities 
in the rotationally-supported disk were lower than in the 
freely-falling envelope. Because the inward velocities in the 
disk are not described by the thin-disk theory, we leave A p 
as an adjustable parameter. The total gas density at any 
point in the model is the sum of the densities in the disk 
and envelope, 



Pgaa — pdisk + Pc 



3.2 The temperature 



(10) 



We assume that the envelope is heated by the star. The en- 
velope temperature is (equation 7.36 of lLamers fc Cassinellil 



T cnv = T t {R t /2r) 



2/(4+p) 



(11) 



where p in our model is an adjustable parameter related to 
the dust opacity and the geometry of the flow. If the den- 
sity structure were spherically symmetric, then p would be 
the exponent in the frequency dependence, v v ', of the Planck 
mean opacity of the dust. In a flattened flow, p can be nega- 
tive if the geometrical dilution of the radiation in the accre- 
tion flow is greater than r~ 2 . This can happen if the dust at 
each radius absorbs and isotropically re-emits the outward 
flowing radiation. Radiation that is emitted perpendicular 
to the disk escapes from the flow. Only the radiation that is 
emitted in the direction along the flattened flow continues 
to heat the dust at larger radii. The result is a decrease in 
the radiation in the disk faster than r~ 2 . 

In the thin-disk theory, the disk is heated by dissipa- 
tion related to t he accret i on ra te. The disk temperature is 
(equation 3.23 of lPringld (|l98ll )) 

I ^R?o )\ V R 



(12) 



We include an adjustable factor, Bt, to allow for additional, 
or possibly less, disk heating. For example, the observations 
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constrain the gas density through the observed optical depth 
of NH3. This means that the density and therefore the ac- 
cretion rate, M, are dependent on the assumed NH3 abun- 
dance which is not well known. The adjustable factor, Bt, 
decouples the disk temperature from the assumed molecu- 
lar abundance. Also the disk temperature may be raised by 
stellar radiation (passive heating). The factor, Bt, allows 
for these effects in an approximate way. 

The gas temperature in the model is the density 
weighted average of the envelope and disk temperatures, 



Teas — 



TdiskPdisk + T cnv p c 
Pdisk + Penv 



(13) 



3.3 The velocity 



The 3 components of the gas velocity in the envelop e are 
given in spherical coordinates by equations 4,5,6 of iKetol 
i|2007l ). but there are error s in ea rlier papers. In par ticular, 
equat ions 5 and 6 of IKetol (|2007T ) a nd equ a tion 8 of lUlrichl 
l| 19761 ) are in correct . Equ ation 8 of lUlrichl (|l976h (same as 
equation 5 of lKetol I2OO3)) does not result in energy conser- 
vation, mv 2 /2 — GM/r, when combined with the other ve- 
locity components. The erro r is small, 1 part in 10 4 . The fo l- 
lowing equations, same as in lMendoza. Canto. Ragal l|2004h . 
are exact. 



v r (r,6) 



ve(r,0) = 



GM» 



1/2 



cos 9 



1/2 



(14) 



1/2 / \ / \ 1/2 

(!M \ ( cos6>o — cos6> \ 1 ^_|_ cos# ^ 



COS tio 



1/2 

, GM, \ sin V I cos ( 
v^r, 9) = I I — — - [ 1 



1/2 



(16) 



r / smu \ cost^o 
The velocity in the disk is simply the Keplerian velocity, 

udisk(i?) = y/GM./R (17) 

where the velocity in the disk is purely azimuthal. We as- 
sume that the radial velocity in the rotationally-supported 
disk is comparatively small. The gas velocity in the model 
is the density weighted average of the envelope and disk ve- 
locities, 



^diskPdisk + fenvp 



(18) 



Pdisk + P env 

where v cnv is given by equations 1141 through ITSI and «disk is 
given by 1171 



3.4 The density singularities in the accretion 
model 

The density in the Ulrich model is singular in the mid-plane 
of the disk at the centrifugal radius and also at the origin. 
These singularities are caused by the convergence of the 
streamlines in the simple mathematical description of the 
flow. This convergence is not expected on physical grounds 
because gas pressure, neglected in the Ulrich model, would 
prevent it. We handle the singularities in two ways. We de- 
fine the computational grid to have an even number of cells 
so that the centers of the middle cells are above and below 
the midplane and around the origin. Second, we smooth the 



density in the radial direction with a Gaussian with a width 
of R d /2. 



3.5 The model parameters 

Based on the analyses of the previous observations cited 
in the introduction, we assume that the disk-envelope is 
viewed edge-on. The model contains 6 adjustable parame- 
ters: 

1) po sets the density of the envelope (equation [1} and the 
mass accretion rate (equation [3} , 

2) p sets the exponent of the power law decrease of the 
temperature in the envelope (equation 1 11[) . 

3) r = v<t,r is the specific angular momentum (equation [4| 
of the envelope flow, 

4) M, is the stellar mass 

5) A p sets the ratio (equation [9j of the disk density to the 
envelope density at Rd, 

6) Bt is factor multiplying the disk temperature (equation 

m. 



4 FITTING THE MODEL TO THE NH 3 DATA 

The disk enters into the model additively, and we can test 
for the presence of a disk by fitting models to the data with 
and without the disk. In the first case, with the disk, we 
adjust all 6 model parameters for both the envelope and the 
disk, and in the second case, without the disk, we adjust 
only the first 4 parameters for the envelope. The fitting is 
done independently in each case, and the 4 envelope param- 
eters are therefore different in the 2 cases. Th e procedure 
for fit ting the data is the same as described in iKeto et al.l 
(2004). We use a fast simulated annealing algorithm to ad- 
just the model parameters to minimize the summed squared 
difference (x 2 ) between the data and the model spectra. 

The model spectra for each particular set of model 
parameter s are generated by our radi ative transfer code 
MOLLIE (|Ketolll990l ; IKeto et alj |2004 ). We assume LTE 
conditions for the NH 3 and CH3CN lines. The LTE ap- 
proximation is appropriate for NH3 because in the absence 
of strong infrared radiation the level populations are ex- 
pected to be mostly in the "metastable" states which are 
the lowest J-state of each K-ladder. Since radiative tran- 
sitions between K-ladders are forbidden, the coupling be- 
tween the metastable states is purely collisional and the 
popul ation approximates a Boltzmann distribution as in 
LTE <|Ho fc Towneslll983l ). CH3CN is also a symmetric top 
and transitions across the K-ladders are similarly forbid- 
den; however, unlike NH 3 , the upper states are easily pop- 
ulated in warm gas (> 100 K). While the justification for 
the LTE approximation is not as strong for CH 3 CN as for 
NH3, we find that CH3CN always traces hot, very dense gas 
where collisional transitions should be important. For C 34 S 
we use the accelerat e d lam bda iteration algorithm (ALI) of 
iRvbicki fc Hummer! (|l99lh to solve for the non-LTE level 
populations. 

In comparing the model to the data, we take into ac- 
count the spatial averaging of the brightness by the width 
of the observing beam. We compute spectra over a grid of 
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positions and smooth the result by convolution with a Gaus- 
sian with the FWHM equivalent to the observing beam of 
each observation, ij2] 

We simultaneously fit 12 NH3 spectra, 4 transitions at 
3 locations. The 3 locations are marked on figure [T] We 
ran twenty thousand trial models for each of the two cases 
with and without the disk. The spectra of the 2 best-fit 
models (with and without the disk) are shown in figures 
[2] and [3] Parameters for both cases are listed in table 1. 
Figure [3] shows the temperature and density in the mid- 
plane of our best fitting models. Figure [5] shows the density 
and velocity of the disk-envelope model on planes parallel 
and perpendicular to the rotation. 

We also have CH 3 CN an d C 34 S data from the observa- 
tions of[c csarom et all i|l999h . The way the radiative trans- 
fer simulation program operates, we cannot use our auto- 
mated search algorithm on more than one molecule at a 
time. The program is recompiled for each molecule. We also 
have not implemented the automated search for CH3CN. 
We could fit the model to the C 34 S data, but this line does 
not have the hyperfine structure of NH3 that is so useful 
in constraining the optical depth and temperature. There- 
fore we opted for a different strategy. We do not use the 
CH3CN or the C 34 S data in the fitting. We use the data 
on these other lines as a check on the model derived from 
the NH3 data. We simulate the CH3CN and C 34 S emission 
using the same two models (with and without the disk) pre- 
viously derived from the NH3 observations for comparison 
with the predicted spectra against the data. The models de- 
rived from the NH3 fitting fix the temperature, densities, 
and velocities, but we still need to assume abundances for 
CH3CN and C 34 S. For CH3CN we assume an abundance of 
6 x 10~ 8 . The assumed abundance of C 34 S is 9.0 x 10" 11 , 
chosen to match the brightness of the (2-1) line. 



4.1 Comparison of the model with the observed 
spectra 

The comparison of figures [2] and [3] shows that the higher 
temperature and density disk and the lower temperature 
and density envelope are both required to fit both the high- 
excitation and low-excitation NH3 spectra. Without the 
disk, there is not enough warm, dense gas to reproduce 
the (4,4) line brightness. Without the large-scale envelope, 
there is not enough cold gas to fit the observed ratios of the 
NH3(1,1) hyperfine lines. 

Further evidence for the disk component comes from the 
CH3CN spectrum. The results again show that the warm, 
high-density, disk component is required to get enough 
brightness in lower frequency (higher velocity) K transi- 
tions. Based on the observed line width, the CH3CN(12-11) 
emission comes from gas very close to the star. The ob- 
served CH3CN(12-11) lines that are not blended with other 
molecular lines, have line widths of 10.0 kms - . The 2 lines 
that appear much broader in the data contain contaminat- 
ing emission from tra nsitions of CH 13 CN(12-11) and HNCO 
(|Cesaroni et al.l 1999). Rotational velocities > 10 kms -1 are 
only found at radii closer than 600 AU (< GM/v 2 ) around 
a 10 Mq star. Figure U shows that the gas in the disk has a 
temperature of greater than 150 K within this radius. Thus 
the CH3CN emission derives mostly from warm gas that is 



very close to the star and is thus a good molecular tracer of 
the inner disk. 

The CH3CN(13-12) spectra is too noisy to help dis- 
criminate between the two cases. The peak signal-to-noise 
ratio is about 6 after smoothing by every other channel. 
Nonetheless, the models are at least consistent with these 
observations. 

The brightness ratio of the C 34 S (5-4) and (2-1) lines 
also suggests the presence of a warm, dense disk (figure [7|). 
The model with the disk reproduces the observed brightness 
ratio, 6:1 for the 2 lines, although line width is less than ob- 
served. The model without the disk is not able to generate 
sufficient brightness in the higher excitation line, and the 
shapes of the line profiles do not match the data. In partic- 
ular, the strong splitting that is seen in the model profiles 
without the disk is not seen in the data. This difference is 
also seen in the NH3 and CH3CN spectra of the 2 mod- 
els, although not as prominently. We usually associate split 
line profile s with spatially unresolved observations of Keple- 
rian disks (jBeckwith fc Sargend[l993T ) whereas here the disk 
produces a triangular profile. What happens here is that the 
disk component is spatially resolved and has a very high den- 
sity. Thus within the beam through the center of the model 
there is a lot of high density gas from the outer part of the 
disk with very low velocity projected along the line-of-sight. 
This gas disk creates the peak in the spectrum around zero 
velocity. In the envelope-only model without the high den- 
sity disk, the model is optically thin in C 34 S. In this case, we 
get the usual split line profiles from the rotation and infall 
in the envelope. 

4.2 Goodness of fit 

How well do the data constrain the model parameters? Some 
appreciation can be gained by plotting \ 2 obtained from all 
the trial models versus each model parameter. Each panel 
shows the fits obtained (ordinate) for each value of a single 
parameter (abscissa) as all the other parameters are varied. 
The lowest value of x 2 (ordinate) at each parameter value 
(abscissa) is the best fit that can be obtained for that param- 
eter value, for any combination of all the other parameters. 
The formal error of a model parameter is proportional to the 
second derivative of x 2 with respect to the model parameter. 
Thus the curvature or width of the lower boundary of the 
collection of points in each figure is a qualitative measure of 
the sensitivity of the model to the parameter. 



5 DISCUSSION 
5.1 The density 

The gas density is constrained by the optical depths of the 
NH3 lines, which can be determined from the brightness 
ratios of the hyperfine lines. The density, po, at Rd is 7.9 x 10 4 
(table 1) assuming an NH3 abundance of 10 -7 . With this 
abundance, the mass of the disk is 2.5 Mq, the mass of the 
envelope is 12.6 Mq, and the accretion rate (equation [3]) , 
asssumed to be the same in the envelope and the disk, is 
M = 7.6 x 10 -5 M yr" 1 . The total mass in the accretion 
flow is mass of the disk and envelope together which is 15.1 
Mq. This mass estimate is the total mass within the model 
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boundary of 26000 AU. ICesaroni et all (| 19991 ) estimates that 
there is between 0.6 and 8.0 M© within a radius of 5000 AU. 

The mass estimate derived from molecular line ob- 
servations is subject to a large uncertainty. The radiative 
transfer modeling determines the column density of NH3 
rather than H2. Therefore, the gas density, the masses of 
the disk and envelope, and the accretion rate depend in- 
versely on the assumed abundance of NH3 which is not well 
known. Estimates from models and observations of simi- 
lar c l ouds range from 10~ 9 to 10~ 6 dHerbst fc Klempererl 
20031; lKetolll990l; lEstalella et al.lll993l ; ICaproni et al.ll2000l ; 



Galvan-Madrid et al.ll2009T ). Thus the abundance of NH. 



may be uncertain by more than an order of magnitude, 
yet a factor of 2 in the masses of the disk and envelope 
is significant in an interpretation of the accretion dynamics. 
Furthermore, the abundance of NH3 could be different in 
the disk and envelope. Estimates of the NH3 abundance are 
generally lowest in colder clouds and highest in warm gas 
around massive stars. Some NH3 may be frozen onto dust 
grains in colder gas and sublimated into the gas phase as 
the temperature rises. Therefore, it is possible that the NH3 
abundance is higher in the warm disk than the cold enve- 
lope. If so, the mass of the envelope could be higher than 
12.6 Mq, assuming a lower abundance, without necessarily 
implying a higher mass for the disk. 



5.2 The disk density factor 

In the model with the flared disk, the density in the mid- 
plane of the disk at the radius of the disk boundary Rd is 5 
times more dense than the smoothed density of the envelope 
at the same point. Considering that the gas density changes 
by 6 orders of magnitude across the model, there is little 
discontinuity between the envelope and disk densities (figure 
[4|. In the model without the disk the smoothed density is 
approximately constant from Rd to the origin (figure [J). 



5.3 The stellar mass 

The stellar mass, 10.7 Mq, is constrained primarily by the 
disk velocities (equation 17 ) required to match the observed 
linewidths (about 10 kms" 1 ) and the gas temperature that 
determines the brightness ratios of the low and high ex- 
citation lines. The gas velocities and therefore line widths 
due to unresolved motions within the beam depend on the 
stellar mass because the velocities are proportional to the 
square root of the mass. The gas temperature depends on 
the mass of the star because the stellar temperature is a 
strong function of the stellar mass and because the accre- 
tion rate depends on the stellar mass. In the disk, heating 
by both the star and by accretion (dissipation) are impor- 
tant (equations 1111 [121 and [T3l ). In the best-fit model, the 
disk temperature would be about 1/3 lower without the pas- 
sive heating from the star, The temperature of the infalling 
envelope is determined entirely by the stellar temperature 
(equation II II) . 



5.4 The angular momentum of the envelope 

The angular momentum of the envelope is determined from 
the widths of the spectral lines and from the VLSR of the 



NH3 spectra at the locations left and right of the center. 
If there were too much rotation then both of the off-center 
spectra would have an incorrect Vlsr; too little, and all 
the NH3 line widths would be too narrow. The disk radius, 
Rd = 6900, is set by the angular momentum in the envelope, 
r, and the stellar mass (equation 2]). The velocity at Rd, 
Vk — 1.2 kms -1 , increases inward as R~ 1 ' 2 . We assume a 
microturbulent broadening of 1 kms -1 , added in quadrature 
to the thermal broadening. Since the observed line widths 
are about 10 kms -1 most of the width of the spectral lines 
comes from rotation and infall in the model. 

The model radius of 6900 AU seems large, and we would 
regard this as an upper limit. First, the Ulrich flow conserves 
angular momentum whereas real accretion flows probably 
involve some braking of the spin-up. If the spin-up were 
slower, then the disk would be found at a smaller radius. 
Second, the high density of the disk is helpful in providing 
optical depth to strengthen the NH3(1,1) outer hyperfines. 
The Ulrich flow has a density singularity in the mid-plane 
which we have handled as desribed in i)3.4l bv a combination 
of gridding and smoothing. Maybe our mid-plane density in 
the envelope is too low, and to compensate, the thin disk is 
bigger. 

5.5 The exponent of the temperature power law 

The modelling suggests that the temperature of the envelope 
falls off very quickly away from the star, as r . This implies 
that the dilution of the radiation is faster than spherical, 
r -2 . The rapid dilution is consistent with the rotationally 
flattened geometry. At each radius, radiation is absorbed 
and re-emitted by dust in the flow. In a flattened flow, much 
of this reprocessed radiation escapes vertically out of the 
flow. In contrast, in a spherical flow, all the reprocessed ra- 
diation continues to interact with the flow at larger radii. 
If the gas temperature decreases outward fast enough that 
most of the envelope is at the minimum gas temperature, 
then a larger value of the exponent (more negative value 
of the parameter p in equation lll[) would have no further 
effect. We assume a minimum temperature of 10 K, a typi- 
cal temperature of molecular gas, We derive an upper limit 
p < -1. 

5.6 The disk temperature multiplier 

The disk temperature multiplier determines the relative im- 
portance of active disk heating, owing to accretion and dissi- 
pation (eouation ll2p . and passive heating from the protostar 
(equation [11). With active disk heating, the temperature in 
the disk is about 50% above that from passive heating alone. 



6 CONCLUSIONS 

This investigation shows that a standard model of a freely- 
falling, rotationally-flattened accretion envelope around a 
rotationally supported disk is able to reproduce the NH3, 
CH3CN, and C 34 S spectral line observations of IRAS20126. 
Both the disk and envelope components are required to fit 
the observed brightness of both the low and high excitation 
lines, and to fit the observed line widths. 

This disk-envelope model was developed for low-mass 
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stars and is quite successful in explaining many of their ob- 
servable characteristics. The success of this model in ex- 
plaining the molecular line observations of the massive star 
IRAS20126 suggests that at least up to 10 M©, the accretion 
processes of massive stars are similar to those of solar mass 
stars. 

There are some differences. Although the mass of the 
disk is uncertain owing to its dependence on the NH3 abun- 
dance, the disk mass is a significant fraction of the mass of 
the star. Furthermore, the extent of the disk is quite large, 
~6900 AU. This suggests that self-gravity in the disk is dy- 
namically important, and therefore, the disk may be unsta- 
ble to local fragmentation and the formation of companion 
stars. 

The physical model in combination with the molecular 
line radiative transfer presented in this paper has a further 
application to other massive star forming regions. The spa- 
tial resolution and sensitivity of the present generation of 
interferometers cannot spatially resolve accretion disks in 
massive protostellar objects at multi-kpc distances. How- 
ever, spectral lines sampling a wide range of densities and 
temperatures can still provide constraints to the physical 
structure of the core and disk. As shown in IRAS20126, the 
spectral lines formed in the higher density and temperature 
regions confirm the presence of an accretion disk. While fu- 
ture observatories such as ALMA will be able to spatially 
resolve the flow on the disk scale, before the science commis- 
sioning of ALMA, this method is a promising technique to 
probe the spatially unresolved regions of the flow and disk 
using data from the current interferometers. 
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Table 1. Model Parameters. 



Parameter 


Symbol 


Disk 


No Disk 


tnv. Density at Hd (cm °) 


PO 


7.9 x 10 


7.9 x 10° 


Temperature power law exp. 


V 


< —1 


0.4 


Angular Momentum (AU kms -1 ) 


r 


8100 


3500 


Stellar Mass (Mq) 


M» 


10.7 


7.3 


Disk Density Ratio 


A P 


5.1 




Disk Temperature factor 


Bj 1 


15.0 




Centrifugal radius (AU) 


R d 


6900 


1900 


Velocity at (kms -1 ) 


Vk 


1.2 


1.8 


Total mass 1 within 0.128 pc (Mq) 




12.6 


10.1 


Disk mass (Mq) 




2.5 




Envelope accretion rate (Mq yr" 1 ) 




7.6 x 10~ 5 


1.0 x 10~ 4 



1 The mass estimates include the envelope and the disk for the case 
with a disk, and the envelope only for no disk. 
Both cases assume an NH3 abundance of 1.0 X 10~ 7 . 
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Figure 1. Images of IRAS20126 in NH 3 (1,1), NH 3 (2,2), and SiO(2-l). The red stars in the NH 3 images show the locations of the NH 3 
spectra in figures [2] and \3\ See text for absolute coordinates. The single star in the SiO image corresponds to the middle star in the other 
images. The angular resolution (FWHM) is indicated by the shaded ellipse, lower-left corner of each panel. The contour interval and 
first contour is 35 mj y kms" 1 for the NH 3 (1,1) and NH 3 (2,2) emission, and 130 mjy kms" 1 for the SiO(2-l). The SiO image is from 
ICesaroni et al.1 dl99STl . 
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Figure 2. Data (red) and model (blue) spectra of NH 3 (1,1) (top row) NH 3 (2,2) (2nd row) NH 3 (3,3) (3rd row) and NH 3 (4,4) (bottom 
row) at 3 positions across the mid-plane of IRAS20126 for the accretion flow including both an infalling envelope and a flared disk. The 
spectra in the middle column are toward the center of the flow, and the the columns to the left and right show spectra on either side of 
the center, 6200 AU to the left, and 5700 au to the right. The 3 locations are shown in figure [T] 
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Figure 3. Spectra in the same format as figure [2] for the accretion flow with an envelope only without a flared disk. 
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Figure 4. The temperature (red) and density (blue) for the model with both an envelope and flared disk (solid line) and for the model 
without a disk (dashed line). The temperatures are about the same in both models, but the disk adds considerable density. The increase 
in density in the envelope-only model (dashed) at 3300 AU is the smoothed density singularity at Rj in the Ulrich (1976) model. 




Figure 5. Density (color) and velocities (arrows) for the model with both an envelope and flared disk. Left: The XZ plane (perpendicular 
to the rotation) and right: the XY plane (the mid- plane of the disk and rotationally- flattened envelope). In the left panel, the color scale 
is logarithmic between 9.2 X 10 4 and 3.7 X 10 6 cm -3 . Contours (white) indicate the density at 0.006, 0.01, and 0.03, of the peak density. 
The longest arrow corresponds to 2.5 kms" 1 In the right panel the color scale is logarithmic between 1.1 X 10 5 and 8.1 X 10 7 cm -3 . 
Contours (white) show the density at 0.03, 0.10, and 0.3 of the peak density. The linear scale in both panels is R^ = 6900 AU. 
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Figure 6. Data (red) and model (blue) spectra of CH3CN(12-11). The left panel shows the spectrum for the model with both an 
accretion envelope and flared disk. The model is the same as shown in figures [2] HI El The right panel is for the model with only an 
accretion envelope and without a disk (figures [3] [Hi ■ In the observed spectrum, the line that is anomalously wide (K=6, 2nd from right) 
is contaminated by HNCO emission. Additionally, contamination from weak CHg 3 CN emission shows up to the red (right) of 4 of the 
lines. The velocity of the K=6 hyperfine line has been set to zero. The observed spectra have been shifted in velocity to match. 




Figure 7. C 34 S(2-1) (blue) and C 34 S(5-4)(red) spectra of the data (solid lines) and the model (dashed lines). Left: Model with both 
the accretion envelope and flared disk. Right: Model with only an accretion envelope and without a flared disk. The brightness of the 
C 34 S(5-4) spectra, both observed and modeled, have been divided by 6. The comparison shows that the model with the disk produces 
the correct brightness ratio. Without the disk, the C 34 S(5-4) line is not bright enough and both lines have the wrong profile shape. The 
zero velocity refers to the model. The observed spectra have been shifted to match. 




-220 -200 -180 
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Figure 8. Spectra of CH3CN(13-12) in the same format as figure [6] except that we have only the 4 lowest K transitions. The K=0 and 
K=l transitions are blended together in the bluest peak. The left panel shows the spectrum for the model with both an accretion envelope 
and flared disk. The right panel is for the model with only an accretion envelope and without a disk (figures [3] [4{ . The signal- to- noise 
ratio of the observed spectra, about 6 at the peak, is not high enough to discriminate between the two models. At least the models are 
consistent with the data. The data are from T.K. Sridharan (private communication). 
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Figure 9. The value of \ 2 as a function of a single parameter. The plots show \ 2 f° r the better fitting disk-envelope models as a function 
of a single parameter. The width of the collection of points in each figure is a qualitative measure of the sensitivity of the models to each 
parameter. Values of the parameters outside the range shown by the lower x values produce very poor fits. The units of the density 
are the log of cm -3 , the units of stellar mass are Mq, and the units of angular momentum arc AU X kms^ 1 divided by 1000. The 
parameters are listed in Table 1. 
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